# import libraries
import pandas as pd
import numpy as np
import seaborn as sns
import plotly.express as px
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime, timedelta
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
from sklearn.model_selection import TimeSeriesSplit
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
import xgboost as xgb
from catboost import CatBoostClassifier
from sklearn.metrics import (
    roc_curve, auc, confusion_matrix, f1_score, precision_score,
    recall_score, accuracy_score, roc_auc_score, classification_report
)
from sklearn.preprocessing import StandardScaler
from matplotlib.dates import YearLocator, DateFormatter

sahm = pd.read_csv('data/SAHM.csv') # Read the CSV file
sahm.columns = ["date", "sahm"]

def create_recession_chart():
    # Create the chart.
    plt.figure(figsize=(12, 6))
    plt.xlabel('Decades')
    plt.grid(True)
    plt.tight_layout()
    # Format the time axis to show only decades.
    decade_fmt = pd.to_datetime([f"{year}-01-01" for year in range(1960, 2030, 10)])
    plt.xticks(decade_fmt, labels=[f"{year}s" for year in range(1960, 2030, 10)])

    # Add 10% blue for recession months.
    recession_periods = {
        (pd.Timestamp('1970-01-01'), pd.Timestamp('1970-11-01')),
        (pd.Timestamp('1973-12-01'), pd.Timestamp('1975-03-01')),
        (pd.Timestamp('1980-02-01'), pd.Timestamp('1980-07-01')),
        (pd.Timestamp('1981-08-01'), pd.Timestamp('1982-11-01')),
        (pd.Timestamp('1990-08-01'), pd.Timestamp('1991-03-01')),
        (pd.Timestamp('2001-04-01'), pd.Timestamp('2001-11-01')),
        (pd.Timestamp('2008-01-01'), pd.Timestamp('2009-06-01')),
        (pd.Timestamp('2020-03-01'), pd.Timestamp('2020-04-01'))
    }
    for start_date, end_date in recession_periods:
        plt.axvspan(start_date, end_date, color='blue', alpha=0.1)

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
sahm['date'] = pd.to_datetime(sahm['date'])
plt.scatter(sahm['date'], sahm['sahm'], s=10, color='gray')  # s = the point size
# Readability settings.
plt.ylabel('SAHM, %')
plt.title('SAHM over time (by decades).')
plt.show()

indpro = pd.read_csv('data/INDPRO.csv') # Read the CSV file
indpro.columns = ["date", "indpro"]

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
indpro['date'] = pd.to_datetime(indpro['date'])
plt.scatter(indpro['date'], indpro['indpro'], s=10, color='green')  # s = the point size
# Readability settings.
plt.ylabel('INDPRO, index')
plt.title('INDPRO over time (by decades).')
plt.show()

sp = pd.read_csv('data/SP500.csv') # Read the CSV file
# Keep only the required columns.
keep_col = ['Date', 'Close']
sp500 = sp[keep_col]
del sp
# Renaming columns.
sp500.columns = ['date', 'sp500']

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
sp500['date'] = pd.to_datetime(sp500['date'])
plt.scatter(sp500['date'], sp500['sp500'], s=10, color='purple')  # s = the point size
# Readability settings.
plt.ylabel('SP500, index')
plt.title('S&P 500 over time (by decades).')
plt.show()

df = pd.read_csv('data/DGS10.csv') # Read the CSV file
df.columns = ["date", "tr10"]

# Replace '.' with NaN and convert column 'tr10' to numeric type.
df['tr10'] = pd.to_numeric(df['tr10'], errors='coerce')

# Group by month and calculate the arithmetic mean for each month.
df['date'] = pd.to_datetime(df['date'])
df.set_index('date', inplace=True)
monthly_avg = df.resample('MS').mean()

# Create a new dataframe with the first number of the month and the arithmetic mean.
tr10 = monthly_avg.reset_index()
del df

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
tr10['date'] = pd.to_datetime(tr10['date'])
plt.scatter(tr10['date'], tr10['tr10'], s=10, color='red')  # s = the point size
# Readability settings.
plt.xlabel('Decades')
plt.ylabel('TR10, %')
plt.title('10 year Treasury Bond over time (by decades).')
plt.show()

data = pd.read_csv('data/T10YFF.csv') # Read the CSV file
data.columns = ["date", "t10yff"]

# Replace '.' with NaN and convert column 't10yff' to numeric type.
data['t10yff'] = pd.to_numeric(data['t10yff'], errors='coerce')

# Group by month and calculate the arithmetic mean for each month.
data['date'] = pd.to_datetime(data['date'])
data.set_index('date', inplace=True)
monthly_avg = data.resample('MS').mean()

# Create a new dataframe with the first number of the month and the arithmetic mean.
t10yff = monthly_avg.reset_index()
del data

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
t10yff['date'] = pd.to_datetime(t10yff['date'])
plt.scatter(t10yff['date'], t10yff['t10yff'], s=10, color='orange')  # s = the point size
# Readability settings.
plt.ylabel('T10YFF, %')
plt.title('Slope of the yield curve over time (by decades).')
plt.show()

unrate = pd.read_csv('data/UNRATE.csv') # Read the CSV file
unrate.columns = ["date", "unrate"]
# => Conversion from % to numbers !!
unrate["unrate"] = unrate["unrate"] / 100

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
unrate['date'] = pd.to_datetime(unrate['date'])
plt.scatter(unrate['date'], unrate['unrate'], s=10, color='brown')  # s = the point size
# Readability settings.
plt.ylabel('UNRATE, %')
plt.title(' Monthly Unemployment Rate over time (by decades).')
plt.show()

pcepi = pd.read_csv('data/PCEPI.csv') # Read the CSV file
pcepi.columns = ["date", "pcepi"]

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
pcepi['date'] = pd.to_datetime(pcepi['date'])
plt.scatter(pcepi['date'], pcepi['pcepi'], s=10, color='blue')  # s = the point size
# Readability settings.
plt.ylabel('PCEPI, index')
plt.title('PCEPI over time (by decades).')
plt.show()

payems = pd.read_csv('data/PAYEMS.csv') # Read the CSV file
payems.columns = ["date", "payems"]

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
payems['date'] = pd.to_datetime(payems['date'])
plt.scatter(payems['date'], payems['payems'], s=10, color='pink')  # s = the point size
# Readability settings.
plt.ylabel('PAYEMS, Thousands of Persons')
plt.title('PAYEMS over time (by decades).')
plt.show()

houst = pd.read_csv('data/HOUST.csv') # Read the CSV file
houst.columns = ["date", "houst"]

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
houst['date'] = pd.to_datetime(houst['date'])
plt.scatter(houst['date'], houst['houst'], s=10, color='magenta')  # s = the point size
# Readability settings.
plt.ylabel('HOUST, Thousands of Units')
plt.title('HOUST over time (by decades).')
plt.show()

recession = pd.read_csv('data/USREC.csv') # Read the CSV file
recession.columns = ["date", "recession"]

px.line(recession, x= recession.date, y='recession')

print("There are a total of {} months of recession in this data.".format(recession["recession"].sum()))

There are a total of 85 months of recession in this data.

# Returns an array where each value is the logarithmic difference between the current
# and previous values ​​in the original array.
'''
    v.shift(1) => shifts the elements in the array down one position,
    with the first element becoming a NaN and the rest being moved forward one position.
    v/v.shift(1) => calculates the ratio of the current value to the previous value in the series.
    np.log(v/v.shift(1)) => calculates the natural logarithm (log to base e) of this ratio.
'''
def log_diff(v): 
    log_diff = np.log(v/v.shift(1))
    return log_diff

# Create a new empty table.
collected_data = pd.DataFrame(columns=[
    'sahm', # Real-time Sahm Rule Recession Indicator
    'indpro', # Industrial Production Index
    'sp500', # S&P 500 Index
    'tr10', # 10 year Treasury Bond
    't10yff', # Slope of the yield curve
    'unrate', # Unemployment Rate
    'pcepi', # Personal Consumption Expenditures
    'payems', # Total Nonfarm Payroll
    'houst', #  New Privately-Owned Housing Units Started
    'recession' # Recession Indicators
])
# Add data to the columns and preprocessing.
collected_data['sahm'] = sahm['sahm']
collected_data['t10yff'] = t10yff['t10yff']
collected_data['unrate'] = unrate['unrate']
collected_data['recession'] = recession['recession']
# Logarithmization is appropriate for series that have exponential growth or high variation.
collected_data['indpro'] = log_diff(indpro['indpro'])
# Logarithmization is appropriate for series that have exponential growth or high variation.
collected_data['sp500'] = log_diff(sp500['sp500'])
# TR10 has a clear trend (increasing or decreasing), differentiation can help.
collected_data['tr10'] = tr10['tr10'].diff()
# PCEPI has a clear trend (increasing or decreasing), differentiation can help.
collected_data['pcepi'] = pcepi['pcepi'].diff()
# PAYEMS has a clear trend (increasing or decreasing), differentiation can help.
collected_data['payems'] = payems['payems'].diff()
# Time series normalization due to large residuals.
# We create an instance of MinMaxScaler.
scaler = MinMaxScaler()
# We apply MinMaxScaler to the 'payems' column.
collected_data['payems'] = scaler.fit_transform(collected_data[['payems']])
# We apply MinMaxScaler to the 'houst' column.
collected_data['houst'] = houst['houst']
collected_data['houst'] = scaler.fit_transform(collected_data[['houst']])

# Drop missing values after preprocessing process.
len_before = len(collected_data)
collected_data.dropna(inplace=True)
len_after = len(collected_data)
print("Dropped {} rows".format(len_before - len_after))

Dropped 1 rows

plt.figure(figsize=(16, 6)) # Determine the size of the figure.
# Define the mask to set the values ​​in the upper triangle to True.
mask = np.triu(np.ones_like(collected_data.corr(method='kendall'), dtype=bool))
heatmap = sns.heatmap(collected_data.corr(), mask=mask, vmin=-1, vmax=1, annot=True, cmap='BrBG')
heatmap.set_title('Triangle Correlation Heatmap', fontdict={'fontsize':18}, pad=16)
plt.show()

# Create list of monthly dates.
dates = pd.date_range(start='02-1962', end='09-2024', freq='MS')
formatted_dates = pd.to_datetime(dates)

# Assuming collected_data is already an existing DataFrame with some data.
# Check the length of collected_data to ensure it matches the length of formatted_dates.
if len(collected_data) == len(formatted_dates):
    collected_data['date'] = formatted_dates
else:
    print("The length of collected_data does not match the length of formatted_dates!")
collected_data.set_index('date', inplace=True)

collected_data # _______ DataFrame (DF) with all data! _______

collected_data.to_csv(r'data/collected_data.csv') # We save the dataframe to a file.

collected_data.describe().T

c1 = collected_data["recession"].value_counts()[0]
c2 = collected_data["recession"].value_counts()[1]
print(collected_data["recession"].value_counts())
print("The minority class only makes up {} % ".format(round(c2/c1 * 100, 2)))

recession
0    667
1     85
Name: count, dtype: int64
The minority class only makes up 12.74 %

# get X and y 
X = collected_data.drop(['recession'], axis=1)
y = collected_data['recession']

# We define the training period.
X_train, y_train = X.loc["1962-02-01":"2012-12-01"], y.loc["1962-02-01":"2012-12-01"]
# We define the test period.
X_test, y_test = X.loc["2013-01-01":], y.loc["2013-01-01":]

# Proper separation of data with stratification.
X_train, X_test, y_train, y_test = train_test_split(
    X, y, 
    test_size=0.2, 
    random_state=42,
    stratify=y  # This is the key change.
)

print("Data form:")
print(f"X_train shape: {X_train.shape}")
print(f"X_test shape: {X_test.shape}")
print(f"y_train shape: {y_train.shape}")
print(f"y_test shape: {y_test.shape}")

print("\nDistribution of classes:")
print("Train class distribution:")
print(pd.Series(y_train).value_counts(normalize=True))
print("\nTest class distribution:")
print(pd.Series(y_test).value_counts(normalize=True))

# Scaling the data.
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Setting the graphics style.
plt.rcParams['figure.figsize'] = (10, 8)
plt.rcParams['axes.grid'] = True

# A function to calculate specificity.
def specificity_score(y_true, y_pred):
    tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()
    return tn / (tn + fp)

# Initialization of models with their hyperparameters.
models = {
    'Logistic Regression': {
        'model': LogisticRegression(
            C=1,
            class_weight={0: 1, 1: 6},
            max_iter=1000,
            penalty='l2',
            solver='saga',
            random_state=42
        ),
        'needs_scaling': True
    },
    'Decision Tree': {
        'model': DecisionTreeClassifier(
            criterion='entropy',
            max_depth=5,
            max_features=None,
            min_samples_leaf=1,
            min_samples_split=2,
            random_state=42
        ),
        'needs_scaling': False
    },
    'Random Forest': {
        'model': RandomForestClassifier(
            bootstrap=False,
            class_weight='balanced',
            max_depth=15,
            max_features='sqrt',
            min_samples_leaf=4,
            min_samples_split=10,
            n_estimators=300,
            random_state=42
        ),
        'needs_scaling': False
    },
    'XGBoost': {
        'model': xgb.XGBClassifier(
            colsample_bytree=0.9,
            gamma=0.1,
            learning_rate=0.05,
            max_depth=6,
            min_child_weight=3,
            n_estimators=100,
            scale_pos_weight=1,
            subsample=1.0,
            random_state=42
        ),
        'needs_scaling': False
    },
    'CatBoost': {
        'model': CatBoostClassifier(
            auto_class_weights='Balanced',
            bagging_temperature=0,
            border_count=32,
            depth=8,
            iterations=200,
            l2_leaf_reg=1,
            learning_rate=0.05,
            random_strength=10,
            random_state=42,
            verbose=False
        ),
        'needs_scaling': False
    }
}

# List to store the results.
results = []
# Creating a figure for the ROC curves.
fig, ax = plt.subplots()

# Training and evaluation of models.
for name, model_info in models.items():
    model = model_info['model']
    needs_scaling = model_info['needs_scaling']
    
    # Choosing the right data according to the scaling need.
    X_train_data = X_train_scaled if needs_scaling else X_train
    X_test_data = X_test_scaled if needs_scaling else X_test
    
    try:
        # Training the model.
        model.fit(X_train_data, y_train)
        
        # Predictions:
        y_pred = model.predict(X_test_data)
        y_pred_proba = model.predict_proba(X_test_data)[:, 1]
        
        # Checking predictions.
        print(f"\nМодел: {name}")
        print("Unique values ​​in predictions:", np.unique(y_pred))
        print("Range of probabilities:", np.min(y_pred_proba), "-", np.max(y_pred_proba))
        
        # Confusion matrix:
        tn, fp, fn, tp = confusion_matrix(y_test, y_pred).ravel()
        print(f"Confusion Matrix: TN={tn}, FP={fp}, FN={fn}, TP={tp}")
        
        # Calculating metrics:
        fpr, tpr, _ = roc_curve(y_test, y_pred_proba)
        roc_auc = auc(fpr, tpr)
        
        # Save metrics.
        results.append({
            'Model': name,
            'AUC': roc_auc,
            'F1-score': f1_score(y_test, y_pred),
            'Precision': precision_score(y_test, y_pred),
            'Recall': recall_score(y_test, y_pred),
            'Specificity': specificity_score(y_test, y_pred)
        })
        
        # Add ROC curve.
        ax.plot(fpr, tpr, lw=2, label=f'{name} (AUC = {roc_auc:.3f})')
        
    except Exception as e:
        print(f"Error at {name}: {str(e)}")

# Finalizing the ROC plot.
ax.plot([0, 1], [0, 1], 'k--', label='Random', lw=2)
ax.set_xlim([0.0, 1.0])
ax.set_ylim([0.0, 1.05])
ax.set_xlabel('False Positive Rate')
ax.set_ylabel('True Positive Rate')
ax.set_title('ROC curves for all models')
ax.legend(loc="lower right")
ax.grid(True)
plt.tight_layout()
plt.show()

# Create a DataFrame with the results.
results_df = pd.DataFrame(results)
results_df = results_df.set_index('Model')
print("\nMetrics for all models:")
print(results_df.round(4).to_string())

Data form:
X_train shape: (601, 9)
X_test shape: (151, 9)
y_train shape: (601,)
y_test shape: (151,)

Distribution of classes:
Train class distribution:
recession
0    0.886855
1    0.113145
Name: proportion, dtype: float64

Test class distribution:
recession
0    0.887417
1    0.112583
Name: proportion, dtype: float64

Модел: Logistic Regression
Unique values ​​in predictions: [0 1]
Range of probabilities: 2.5377559223161707e-07 - 0.999983908938828
Confusion Matrix: TN=122, FP=12, FN=0, TP=17

Модел: Decision Tree
Unique values ​​in predictions: [0 1]
Range of probabilities: 0.0 - 1.0
Confusion Matrix: TN=130, FP=4, FN=2, TP=15

Модел: Random Forest
Unique values ​​in predictions: [0 1]
Range of probabilities: 0.0 - 0.9962633130021159
Confusion Matrix: TN=130, FP=4, FN=3, TP=14

Модел: XGBoost
Unique values ​​in predictions: [0 1]
Range of probabilities: 0.002605706 - 0.9757932
Confusion Matrix: TN=132, FP=2, FN=2, TP=15

Модел: CatBoost
Unique values ​​in predictions: [0 1]
Range of probabilities: 7.03488963512961e-05 - 0.9993221405581604
Confusion Matrix: TN=130, FP=4, FN=1, TP=16

Metrics for all models:
                        AUC  F1-score  Precision  Recall  Specificity
Model                                                                
Logistic Regression  0.9794    0.7391     0.5862  1.0000       0.9104
Decision Tree        0.9306    0.8333     0.7895  0.8824       0.9701
Random Forest        0.9820    0.8000     0.7778  0.8235       0.9701
XGBoost              0.9868    0.8824     0.8824  0.8824       0.9851
CatBoost             0.9947    0.8649     0.8000  0.9412       0.9701

# We create a DataFrame with predictions using the X_train index.
predictions_df = pd.DataFrame(index=X_train.index)
# We add the predictions from each model.
for name, model_info in models.items():
    model = model_info['model']
    needs_scaling = model_info['needs_scaling']
    
    X_data = X_train_scaled if needs_scaling else X_train
    predictions = model.predict_proba(X_data)[:, 1]
    
    predictions_df[name] = predictions

# We sort the index to ensure it is monotonic.
predictions_df = predictions_df.sort_index()

# We find the closest date to the beginning and end of the desired period.
start_date = predictions_df.index[predictions_df.index >= '1990-01-01'].min()
end_date = predictions_df.index[predictions_df.index <= '2012-12-31'].max()

# We filter the data only for the period 1990-2012.
predictions_df = predictions_df.loc[start_date:end_date]

# We create the graphic.
plt.figure(figsize=(15, 10))
# Colors for every model.
colors = ['blue', 'green', 'red', 'purple', 'orange']

# We visualize the predictions for each model.
for (name, color) in zip(models.keys(), colors):
    plt.plot(predictions_df.index, predictions_df[name], color=color, label=name)

# We add vertical columns for the specific periods.
periods = [
    ('1990-08-01', '1991-03-01'),
    ('2001-04-01', '2001-11-01'),
    ('2008-01-01', '2009-06-01')
]

for start, end in periods:
    plt.axvspan(pd.to_datetime(start), pd.to_datetime(end), color='blue', alpha=0.1)

plt.title(f'RecessionProbability ({start_date.strftime("%Y-%m")} until {end_date.strftime("%Y-%m")})')
plt.xlabel('Year')
plt.ylabel('Probability of recession')
plt.legend()
plt.grid(True)

# Format x-axis to show only years.
plt.gca().xaxis.set_major_locator(YearLocator())
plt.gca().xaxis.set_major_formatter(DateFormatter('%Y'))
plt.tight_layout()

plt.show()

collected_data = pd.read_csv('data/collected_data.csv') # Read the CSV file
data = collected_data.copy()
# Setting column 'date' as index.
data = data.set_index('date', drop=True)
# get X and y 
X = data.drop(['recession'], axis=1)
y = data['recession']

# We define the training period.
X_train, y_train = X.loc["1962-02-01":"2012-12-01"], y.loc["1962-02-01":"2012-12-01"]
# We define the test period.
X_test, y_test = X.loc["2013-01-01":], y.loc["2013-01-01":]

# Building and training the model.
cb_model = CatBoostClassifier(
    auto_class_weights='Balanced',
    iterations=200,
    depth=8,
    learning_rate=0.05,
    random_state=42,
    verbose=False
)

# Training the model.
cb_model.fit(X_train, y_train)

# Using the prediction model.
predictions = cb_model.predict(X_test)

# Recession probabilities (ie, how likely a recession is 6 months from now).
probabilities = cb_model.predict_proba(X_test)[:, 1]  # Probability of Class "1" (Recession)

# Model evaluation:
accuracy = accuracy_score(y_test, predictions)
precision = precision_score(y_test, predictions)
recall = recall_score(y_test, predictions)
f1 = f1_score(y_test, predictions)
roc_auc = roc_auc_score(y_test, probabilities)

print(f"Accuracy: {accuracy:.4f}")
print(f"Precision: {precision:.4f}")
print(f"Recall: {recall:.4f}")
print(f"F1 Score: {f1:.4f}")
print(f"ROC AUC: {roc_auc:.4f}")

print("\nConfusion Matrix:")
print(confusion_matrix(y_test, predictions))

print("\nClassification Report:")
print(classification_report(y_test, predictions))

# ROC curve visualization.
fpr, tpr, thresholds = roc_curve(y_test, probabilities)
roc_auc = auc(fpr, tpr)

plt.figure()
plt.plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (area = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic')
plt.legend(loc="lower right")
plt.show()

Accuracy: 0.9645
Precision: 0.2857
Recall: 1.0000
F1 Score: 0.4444
ROC AUC: 0.9928

Confusion Matrix:
[[134   5]
 [  0   2]]

Classification Report:
              precision    recall  f1-score   support

           0       1.00      0.96      0.98       139
           1       0.29      1.00      0.44         2

    accuracy                           0.96       141
   macro avg       0.64      0.98      0.71       141
weighted avg       0.99      0.96      0.97       141

# We use the latest available data (for example, the last 6 months).
# This will be our current "input" for the model to make a forecast 6 months ahead.
latest_data = X.tail(6)  # We take the last 6 months of available data.

# We predict the probability of a recession 6 months ahead.
future_probabilities = cb_model.predict_proba(latest_data)[:, 1]  # Probability of class "1" (recession).

# We create the correct future dates starting from October 2024.
start_date = pd.Timestamp('2024-10-01')
future_dates = pd.date_range(start=start_date, periods=6, freq='MS')  # MS means start of month.

# Visualization of the likelihood of a recession in the next 6 months.
plt.figure(figsize=(10, 6))

# We are adding color zones for the different risk levels.
plt.axhspan(0.7, 1.0, color='red', alpha=0.1, label='High risk (>70%)')
plt.axhspan(0.3, 0.7, color='yellow', alpha=0.1, label='Medium risk (30-70%)')
plt.axhspan(0, 0.3, color='green', alpha=0.1, label='Low risk (<30%)')

plt.plot(future_dates, future_probabilities, marker='o', linestyle='-', color='b', label='Probability of recession')

# Format x-axis dates.
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))

plt.xlabel('Date')
plt.ylabel('Probability of recession')
plt.title('Probability of a recession in the next 6 months.')
plt.xticks(future_dates, rotation=45)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# Probability of a recession in the next months.
future_prob_df = pd.DataFrame({
    'date': future_dates,
    'probability_of_recession': future_probabilities
})
future_prob_df.set_index('date', inplace=True)

future_prob_df

# Probability of recession => as arithmetic mean of probabilities for all months.
average_probability = (future_prob_df['probability_of_recession'].mean() * 100).round(2)
formatted_probability = f"{average_probability} %"
print(">>> ------ >> ----- > ----")
print(f"Average probability of a recession over the entire period: {formatted_probability}")

>>> ------ >> ----- > ----
Average probability of a recession over the entire period: 71.23 %

new_data = pd.read_csv('data/combined_data.csv') # We read the CSV file with the newly generated data.
# Setting column 'Unnamed: 0' as index.
new_data = new_data.rename(columns={'Unnamed: 0': 'date'})
# We make the 'date' column an ​​index.
new_data = new_data.set_index('date')

new_data

# get X and y 
X = new_data.drop(['recession'], axis=1)
y = new_data['recession']

# We define the training period.
X_train, y_train = X.loc["1962-02-01":"2012-12-01"], y.loc["1962-02-01":"2012-12-01"]
# We define the test period.
X_test, y_test = X.loc["2013-01-01":], y.loc["2013-01-01":]

# Building and training the model.
new_cb_model = CatBoostClassifier(
    auto_class_weights='Balanced',
    iterations=200,
    depth=8,
    learning_rate=0.05,
    random_state=42,
    verbose=False
)

# Training the model.
new_cb_model.fit(X_train, y_train)

<catboost.core.CatBoostClassifier at 0x25c35bc5730>

# We make predictions for the entire test period.
all_probabilities = new_cb_model.predict_proba(X_test)[:, 1]

# We create a DataFrame with all predictions.
full_results_df = pd.DataFrame({
    'Date': X_test.index,
    'Recession Probability': all_probabilities
})

# We only take the last 6 months for preview.
results_df = full_results_df.tail(6)

# We create the visualization.
plt.figure(figsize=(12, 6))
sns.set_style("whitegrid")
plt.plot(results_df['Date'], results_df['Recession Probability'], 
         marker='o', linewidth=2, markersize=8, color='blue')

# We are adding color zones for the different risk levels.
plt.axhspan(0.7, 1.0, color='red', alpha=0.1, label='High risk (>70%)')
plt.axhspan(0.3, 0.7, color='yellow', alpha=0.1, label='Medium risk (30-70%)')
plt.axhspan(0, 0.3, color='green', alpha=0.1, label='Low risk (<30%)')

# Chart formatting.
plt.title('Recession Probability Forecast\nSeptember 2024 - February 2025', 
          fontsize=14, pad=20)
plt.xlabel('Date', fontsize=12)
plt.ylabel('Probability', fontsize=12)

# Format the y-axis as a percentage.
plt.gca().yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: '{:.0%}'.format(y)))
# Date rotation.
plt.xticks(rotation=45)
plt.legend(loc='upper right', title='Levels of risk')

# We add annotations with the exact values.
for idx, row in results_df.iterrows():
    plt.annotate(f'{row["Recession Probability"]:.1%}',
                (idx, row["Recession Probability"]),
                textcoords="offset points",
                xytext=(0,10),
                ha='center',
                fontweight='bold')

# We set the borders of the graph.
plt.ylim(-0.05, 1.05)
# Layout optimization.
plt.tight_layout()
plt.show()

results_df

# Probability of recession => as arithmetic mean of probabilities for all months.
avg_probability = (results_df['Recession Probability'].mean() * 100).round(2)
f_probability = f"{avg_probability} %"
print(">>> ------ >> ----- > ----")
print(f"Average probability of a recession over the entire period: {f_probability}")

>>> ------ >> ----- > ----
Average probability of a recession over the entire period: 35.05 %

def rolling_window_prediction(df, features, target, model, window_size=24, prediction_horizon=6):
    """
    Rolling window approach for recession prediction.
    
    Parameters:
    df : pandas DataFrame
        Input data with features and target variable.
    features : list
        List of feature names.
    target : str
        Target variable (e.g., 'recession').
    model : sklearn-like model
        Model to be trained and used for predictions.
    window_size : int
        Size of the rolling window for training (in months).
    prediction_horizon : int
        Number of months to predict ahead.
    
    Returns:
    pandas DataFrame
        Results with dates, predicted probabilities, and risk levels.
    """
    predictions = []
    probabilities = []
    dates = []
    
    for i in range(len(df) - window_size - prediction_horizon + 1):
        # Define training and future windows
        train_window = df.iloc[i:i + window_size]
        future_window = df.iloc[i + window_size:i + window_size + prediction_horizon]
        
        # Train model on the rolling window
        X_train = train_window[features]
        y_train = train_window[target]
        
        # Skip windows where the target has only one unique value
        if y_train.nunique() <= 1:
            continue
        
        model.fit(X_train, y_train)
        
        # Predict on future window
        X_future = future_window[features]
        pred = model.predict(X_future)
        prob = model.predict_proba(X_future)[:, 1]
        
        predictions.extend(pred)
        probabilities.extend(prob)
        dates.extend(future_window.index)
    
    # Create a results DataFrame
    results = pd.DataFrame({
        'date': dates,
        'predicted_recession': predictions,
        'recession_probability': probabilities
    }).set_index('date')
    
    # Add risk levels
    results['risk_level'] = pd.cut(results['recession_probability'], 
                                   bins=[0, 0.3, 0.6, 1], 
                                   labels=['Low', 'Medium', 'High'])
    return results

# Load and prepare data
data = pd.read_csv('data/collected_data.csv').set_index('date')
data.index = pd.to_datetime(data.index)

# Define features and target
features = [col for col in data.columns if col != 'recession']
target = 'recession'

# Create the model
cboost_model = CatBoostClassifier(
    auto_class_weights='Balanced',
    iterations=200,
    depth=8,
    learning_rate=0.05,
    random_state=42,
    verbose=False
)

# Use rolling window prediction
results = rolling_window_prediction(data, features, target, cboost_model, window_size=24, prediction_horizon=6)

# Visualize results
fig, ax = plt.subplots(figsize=(15, 8))
results['recession_probability'].plot(ax=ax, label='Recession Probability', color='red')
results[results['predicted_recession'] == 1].index.map(
    lambda x: ax.axvline(x=x, color='blue', linestyle='--', alpha=0.1, label='Predicted Recession' if x == results.index[0] else "")
)
ax.set_title('Recession Predictions (Rolling Window Approach)', fontsize=14)
ax.set_xlabel('Date', fontsize=12)
ax.set_ylabel('Probability', fontsize=12)
ax.legend()
plt.tight_layout()
plt.show()

# Parameters for the new DataFrame.
dates = pd.date_range(start='2024-10-01', end='2025-03-01', freq='MS')  # Monthly index

# Calculating the probability of a recession by taking only the last 6 values.
recession_probabilities = (
    future_prob_df['probability_of_recession'].iloc[-6:].values * 0.7 * 100 + 
    results_df['Recession Probability'].iloc[-6:].values * 0.3 * 100
)

# Create the new DataFrame.
new_df = pd.DataFrame(data=recession_probabilities, index=dates, columns=['Recession Probability'])

# Create a graph of the probability of a recession.
plt.figure(figsize=(10, 6))
plt.plot(new_df.index, new_df['Recession Probability'])
plt.title('Probability of recession from September 2024 to February 2025 - ensemble of the two models.')

# We are adding color zones for the different risk levels.
plt.axhspan(70, 100, color='red', alpha=0.1, label='High risk (>70%)')
plt.axhspan(30, 70, color='yellow', alpha=0.1, label='Medium risk (30-70%)')
plt.axhspan(0, 30, color='green', alpha=0.1, label='Low risk (<30%)')

plt.xlabel('Date')
plt.ylabel('Probability of recession. [ % ]')
plt.grid()
plt.show()

new_df

# Probability of recession => as arithmetic mean of probabilities for all months.
avg_probability = (new_df['Recession Probability'].mean()).round(2)
f_probability = f"{avg_probability} %"
print(">>> ------ >> ----- > ----")
print("Average recession probability of the ensemble of the three methods over the entire period: ", f_probability)

>>> ------ >> ----- > ----
Average recession probability of the ensemble of the three methods over the entire period:  60.37 %

create_recession_chart()
# We make sure that the 'date' column is of type datetime.
indpro['date'] = pd.to_datetime(indpro['date'])
plt.scatter(indpro['date'], indpro['indpro'], s=10, color='green')  # s = the point size
# Readability settings.
plt.ylabel('INDPRO, index')
plt.title('INDPRO over time (by decades).')
# Add a horizontal line.
plt.axhline(y=103, color='r', linestyle='--', xmin = 0.05, xmax = 0.68, linewidth=0.5)
plt.axhline(y=103, color='r', linestyle='--', xmin = 0.70, xmax = 0.95, linewidth=2)
plt.show()

Type (bucket)	Base indicator	Index	Modification
Indicator (FED)	Real-time Sahm Rule Recession Indicator	SAHMREALTIME	12 month change
All Markets	Personal Consumption Expenditures	PCEPI	12 month change
Commodity market	Industrial Production Index	INDPRO	12 month change
Stock market	S&P 500 Index	SP500	12 month change
Bond market	10 year Treasury Bond	TR10	12 month change
Bond market	Slope of the yield curve	T10YFF	12 month change
Labor market	Total Nonfarm Payroll	PAYEMS	12 month change
Housing market	New Privately-Owned Housing Units Started	HOUST	12 month change
Unemployment	Monthly Unemployment Rate	UNRATE	12 month change

	sahm	indpro	sp500	tr10	t10yff	unrate	pcepi	payems	houst	recession
date
1962-02-01	-0.17	0.016229	0.016139	-0.043737	1.650556	0.055	0.042	0.827913	0.396825	0
1962-03-01	-0.17	0.005350	-0.005878	-0.108990	1.078182	0.056	0.021	0.819544	0.478671	0
1962-04-01	-0.10	0.002130	-0.063973	-0.087455	1.043000	0.056	0.018	0.829109	0.518849	0
1962-05-01	-0.07	-0.001066	-0.089914	0.030636	1.441818	0.055	0.010	0.817113	0.498512	0
1962-06-01	0.00	-0.002132	-0.085381	0.035411	1.206667	0.055	0.010	0.816714	0.459325	0
...	...	...	...	...	...	...	...	...	...	...
2024-05-01	0.37	0.006954	0.046904	-0.056818	-0.847727	0.040	-0.010	0.824685	0.415179	0
2024-06-01	0.43	0.001500	0.034082	-0.177010	-1.024737	0.041	0.145	0.820780	0.422123	0
2024-07-01	0.53	-0.006214	0.011258	-0.056627	-1.081364	0.043	0.202	0.821816	0.388889	0
2024-08-01	0.57	0.003373	0.022578	-0.377727	-1.459091	0.042	0.142	0.822414	0.437996	0
2024-09-01	0.50	-0.002832	0.019996	-0.147409	-1.406500	0.041	0.217	0.826200	0.434524	0

	count	mean	std	min	25%	50%	75%	max
sahm	752.0	0.412527	0.984315	-0.370000	-0.030000	0.030000	0.370000	9.500000
indpro	752.0	0.001897	0.009524	-0.142045	-0.001726	0.002446	0.006304	0.063771
sp500	752.0	0.005887	0.043518	-0.245428	-0.017919	0.009800	0.034384	0.151043
tr10	752.0	-0.000478	0.277827	-1.755317	-0.145981	-0.002482	0.149554	1.612464
t10yff	752.0	0.982359	1.619839	-6.291905	0.118571	1.157828	2.181648	3.803636
unrate	752.0	0.058871	0.017200	0.034000	0.046000	0.056000	0.070000	0.148000
pcepi	752.0	0.143822	0.156321	-1.058000	0.057750	0.131000	0.209000	1.089000
payems	752.0	0.821600	0.032212	0.000000	0.818827	0.823410	0.826967	1.000000
houst	752.0	0.475000	0.192394	0.000000	0.358879	0.489087	0.584325	1.000000
recession	752.0	0.113032	0.316843	0.000000	0.000000	0.000000	0.000000	1.000000

Type of periods	Start	Stop	Recessions
Train	02.1962	12.2012	Yes, 7
Test	01.2013	09.2024	Yes, 1

Model	F1-score	Precision	Recall	Accuracy
Logistic Regression	0.8024	0.7976	0.8072	0.9460
Decision Tree	0.8075	0.8333	0.7831	0.9493
Random Forest	0.8842	0.8926	0.8787	0.9689
XGBoost	0.8557	0.9267	0.8088	0.9657
CatBoost	0.8675	0.8575	0.9051	0.9624
SVM	0.7677	0.7810	0.7978	0.9330

Poject: Analysis of the US Economy - part 2¶

Author: @engine, Date: November 2024, Machine Learning project¶

Abstract:¶

1. Defining the problem and main assumptions¶

2. Preparing the Data¶

2.1 Criteria according to which the data for the present study were collected.¶

2.2 Feature Selection.¶

2.3 Data collection.¶

3. Modeling¶

3.1 Class Imbalance¶

3.2 Modeling¶

3.3 Testing the models¶

4. Final forecast¶

4.1 First way to forecast.¶

4.2 Second way of prediction.¶

4.3 Third Way of Forecasting.¶

4.4 Final lines.¶

5. Conclusions¶

Resources:¶

Base indicator	Index	Conclusion
Real-time Sahm Rule Recession Indicator	SAHM	stationary
Industrial Production Index	INDPRO	stationary
S&P 500 Index	SP500	stationary
10 year Treasury Bond	TR10	stationary
Slope of the yield curve	T10YFF	stationary
Monthly Unemployment Rate	UNRATE	stationary
Personal Consumption Expenditures	PCEPI	stationary
Total Nonfarm Payroll	PAYEMS	stationary
New Privately-Owned Housing Units Started	HOUST	stationary

	probability_of_recession
date
2024-10-01	0.873895
2024-11-01	0.059752
2024-12-01	0.682443
2025-01-01	0.979363
2025-02-01	0.717093
2025-03-01	0.961044

	sahm	indpro	sp500	tr10	t10yff	unrate	pcepi	payems	houst	recession
date
1962-02-01	-0.170000	0.016229	0.016139	-0.043737	1.650556	0.055000	0.042000	0.827913	0.396825	0
1962-03-01	-0.170000	0.005350	-0.005878	-0.108990	1.078182	0.056000	0.021000	0.819544	0.478671	0
1962-04-01	-0.100000	0.002130	-0.063973	-0.087455	1.043000	0.056000	0.018000	0.829109	0.518849	0
1962-05-01	-0.070000	-0.001066	-0.089914	0.030636	1.441818	0.055000	0.010000	0.817113	0.498512	0
1962-06-01	0.000000	-0.002132	-0.085381	0.035411	1.206667	0.055000	0.010000	0.816714	0.459325	0
...	...	...	...	...	...	...	...	...	...	...
2024-11-01	0.371263	0.001896	-0.000628	0.000468	-1.054334	0.042184	0.199241	0.820944	0.434200	0
2024-12-01	0.365931	0.001896	0.012222	-0.000439	-0.907669	0.042725	0.196893	0.821594	0.430898	0
2025-01-01	0.369188	0.001896	-0.000270	-0.000439	-0.771676	0.043291	0.194642	0.821594	0.435397	0
2025-02-01	0.373684	0.001896	0.011875	-0.000439	-0.645579	0.043803	0.192485	0.821594	0.432419	0
2025-03-01	0.377939	0.001896	0.000068	-0.000439	-0.528658	0.044322	0.190417	0.821594	0.436556	0

	Date	Recession Probability
141	2024-10-01	0.862716
142	2024-11-01	0.476150
143	2024-12-01	0.187531
144	2025-01-01	0.219110
145	2025-02-01	0.168936
146	2025-03-01	0.188394

	Recession Probability
2024-10-01	87.054083
2024-11-01	18.467138
2024-12-01	53.396982
2025-01-01	75.128673
2025-02-01	55.264570
2025-03-01	72.924928